%Solve model allowing for alpha\not=0

%Create grid for q; solve model given q
%Procedure analogous to OptContract.m
qSpace=linspace(0,qBar,qGrid);

for z=1:qGrid
    
    q=qSpace(z);
    
    %Calculate second best outcome, i.e., absent moral hazard over
    %screening
    SBAlpha;
    dFStart=0;
    
    %Set boundary condition for F'(V^B)
    if (aB<0.0001)

        dFStart=(FB*(1+alpha*q)^2-phi*(gamma-r))/phi;

    end
    
    right=1;
    
    if (kappa*q<VB)
        
        right=0;
        
    end
    %Solve HJB equation
    sol = ode15s(@ffunAlpha,[VB max(kappa*q,VB+0.001)+0.001*(kappa*q==VB)],[FB dFStart+tolerance-2*tolerance*right],options,alpha, r, gamma ,Lambda, phi,q,aBar,delta);

    
    %Evaluate HJB equation at V_0
    F=sol.y(1,end);
    dF=sol.y(2,end);

    %Determine initial Monitoring
    aEnd=(F-dF*(sol.x(end)+phi)-phi*(gamma-r))/phi;
    
    
    
    [m,ind]=min(abs(sol.x-kappa*q));
    V0=sol.x(ind);
    F0=sol.y(1,ind);
    
    %Save output in vector

    %Calculate F_{0^-}
    vF0(z)=F0-0.5*kappa*q^2;
    vAB(z)=aB;
    vAEnd(z)=aEnd;
    vV0(z)=V0;
    %Stop iteration when degenerate outcome is achieved
    if (vF0(z)<-10)

        vF0(z+1:qGrid)=zeros(qGrid-z,1);

        break;

    end
    
    
end

%Determine optimal q
[m,ind]=max(vF0);
q=qSpace(ind);
a0=vAEnd(ind);
V0=vV0(ind);

%Solve model under optimal q

    SBAlpha;
    dFStart=0;

    if (aB<0.0001)

        dFStart=(FB*(1+alpha*q)^2-phi*(gamma-r))/phi;

    end
    
    right=1;
    
    if (kappa*q<VB)
        
        right=0;
        
    end
    
    sol = ode15s(@ffunAlpha,[VB max(kappa*q,VB+0.001)+0.001*(kappa*q==VB)],[FB dFStart+tolerance-2*tolerance*right],options,alpha, r, gamma ,Lambda, phi,q,aBar,delta);
    %Calculate state-dependent model quantities
    for i=1:length(sol.x)

            
            V=sol.x(i);
            F1=sol.y(1,i);
            dF1=sol.y(2,i);
            a=(1+alpha*q)*(F1-dF1*V)-phi*(gamma-r)/(1+alpha*q)-dF1*phi/(1+alpha*q);
            a=a/max(0.001,phi+2*dF1*alpha*phi/(1+alpha*q));
        
            if (a<0)
            
                a=0;
          
            
            end
        
            if (a>aBar)
            
                a=aBar;
           
            
            end

            W=phi*a/(1+alpha*q);


            if (W>F1)

                W=F1;
                a=(1+alpha*q)*W/phi;

            end

            lambda=Lambda-a-q-alpha*a*q;
            dV=(gamma+lambda+delta)*V-W*(1+alpha*a);
            vA(i)=a;
            vW(i)=W;



    end
   %Determine initial monitoring effort
    a0=vA(end);
